#ifndef CLASSICALRK_HPP
#define CLASSICALRK_HPP

#include <iostream>
#include <vector>
#include <math.h>
#include <fstream>
#include "Function.hpp"

using namespace std;

class ClassicalRK : public TimeIntegrator
{
private:
    int n_steps;
    double T, mu;
    vector<double> value;
    vector<vector<double>> result;
    Norm norm;
public:
    ClassicalRK() {}
    ClassicalRK(int n_, double time, vector<double> Val) : n_steps(n_), T(time), value(Val) {}
    ~ClassicalRK() {}
    void solver();
    void error();
};

void ClassicalRK::solver()
{
    ofstream outfile;
    outfile.open("output.m");
    outfile << "points = [";
    double k = 1.0 * T / n_steps;
    vector<double> U = value;
    outfile << U << endl;
    result.push_back(U);
    for (int i = 0; i < n_steps; i++)
    {
        auto y1 = f(U, i * k);
        auto y2 = f(U + (0.5 * k) * y1, i * k + 0.5 * k);
        auto y3 = f(U + (0.5 * k) * y2, i * k + 0.5 * k);
        auto y4 = f(U + k * y3, i * k + k);
        U = U + (1.0 * k / 6) * (y1 + 2.0 * y2 + 2.0 * y3 + y4);
        outfile << U << endl;
        result.push_back(U);
    }
    outfile << "];" << endl;
    outfile << "scatter(points(:,1), points(:,2),\".\");" << endl;
    outfile.close();
}

void ClassicalRK::error()
{
    // cout << result[n_steps] << endl;
    // cout << result[0] << endl;
    double error = norm.Inftynorm(result[n_steps] - result[0]);
    cout << error << endl;
}

TimeIntegrator* createClassicalRK(int n_, int p_, double T_, vector<double> value_) { return new ClassicalRK(n_, T_, value_); }
#endif